
##################################################
####### Appendix Figures B1, B2, D1 and E1 #######
##################################################

library(foreign)
library(lme4)
library(scales)
library(ggplot2)
library(memisc)
library(Zelig)
library(MASS)
library(locfit)
library(dplyr)

###############################################################
# Section B: Government Fragmentation - Cross-National Figures 
###############################################################

#####################
##### Figure B1 #####
#####################

lgraph = function(X,Y,CONDITION, main="", xlab="X", ylab="Y", ylim=NULL, deg=2, ...){
  fit <- locfit(Y[CONDITION]~lp(X[CONDITION], deg=deg, ...)); crit(fit) <- crit(fit,cov=0.95); 
  plot(fit,band="local", main=main, xlab=xlab, ylab=ylab, ylim=ylim)
}

AdminChange=read.dta("~services_admin_tscs.dta", convert.factors=FALSE)
str(AdminChange)

row_to_keep = which(AdminChange$wdi_pop>500000)
AdminChangeP = AdminChange[row_to_keep,]

RegionMeans <- AdminChangeP %>% 
  group_by(year, region_x) %>%
  summarise(adminpcM = mean(adminpc), 
            adminM=mean(admin)) 

AfricaMean= filter(RegionMeans, region_x=="Sub-Saharan Africa")
attach(AfricaMean)

# Regional Trajectory -- admin units
#pdf(file="~AfricaAdminChangeOT.pdf", width=12, height=10) 
deg=3
lgraph(year, adminM, TRUE, main="Sub-Saharan Africa", xlab="", ylab="Top-tier local governments", deg=deg)

#####################
##### Figure B2 #####
#####################

AfricaAdminChange= AdminChangeP  %>% 
  filter(region_x=="Sub-Saharan Africa" & cname!="Uganda" & cname!="Cape Verde")  %>% 
  select(year, admin, adminpc, cname) %>% 
  arrange(cname, year)

#pdf(file="~AfricaAdminChangeByCountry.pdf", width=12, height=10) 
ggplot(AfricaAdminChange, aes(y=admin, x=year)) +geom_line(colour="red") +
  facet_wrap(~cname)+ theme_bw() +
  scale_x_continuous(name="", breaks = seq(1960, 2010, 20))  

################################################################
# Section D: Robustness checks: Country-level Analysis 
################################################################

#####################
##### Figure D1 #####
#####################

#models_all.dta produced by Appendix_Cross_National_Replication_Stata.do

data <- read.dta("~models_all.dta")

data <- subset(data,data$parm=="ladminpc_l5")
data$model <- seq(1:46)

#pdf(file="~/Dropbox/AdministrativeChanges/draft/appendix/jackknife.pdf",width=10)
ggplot(data = data, aes(x = factor(model), y = estimate)) +geom_point(shape=21, size=3, fill="black")+
  ggtitle("") +
  geom_errorbar(aes(ymin=min95, ymax=max95), width=0.05, size=0.5) +
  scale_shape_manual(values=c(22,22))  + scale_y_continuous(name="Regression Coefficient",breaks = round(seq(min(0), max(0.8), by = 0.05),3), limits = c(0, 0.6)) + 
  scale_x_discrete(name="Model") +  scale_fill_manual(values=c("black","black")) + coord_flip() +
  geom_hline(yintercept=0, linetype="longdash", size=1) +  theme_bw(base_size = 12, base_family = "Helvetica") +
  theme(legend.position=c(0.9, 0.5),
        panel.border = element_blank(), panel.grid.major.y = element_line(),
        panel.grid.major.x = element_blank(), panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"),plot.title = element_text(size = 16, face = "bold", vjust = 1.5),
        axis.title.y = element_text(colour="grey20",size=14,face="bold"),
        axis.text.x = element_text(colour="grey20",size=13, face="bold"),
        axis.text.y = element_text(colour="grey20",size=12,face="bold"),  
        axis.title.x = element_text(colour="grey20",size=14,face="bold"), 
        axis.ticks.y = element_blank())


################################################################
# Section E: Section E: Sub-national analysis 
################################################################

#####################
##### Figure E1 #####
#####################

# real data
setwd("~")
admin <- read.dta("~services_admin_tscs.dta")
admin <- admin[admin$region_x=="Sub-Saharan Africa",]

### only looking at the effect until 10 units pc
draws <- mvrnorm(n=1000, means, var)
max <- 10
adminvec <- seq(min(admin$adminpc,na.rm=T),max, by=((max)-min(admin$adminpc,na.rm=T))/(500-1))

meanvec <- lowervec <- highervec <- NULL
for(i in 1:length(adminvec)){
  out <- draws[,43] + varmeans[1]*draws[,1]+ varmeans[2]*draws[,2]+varmeans[3]*draws[,3]+varmeans[4]*draws[,4]+varmeans[5]*draws[,5]+varmeans[6]*draws[,6]+varmeans[7]*draws[,7]+varmeans[8]*draws[,8]+adminvec[i]*draws[,9]+ (adminvec[i]^2)*draws[,10]+varmeans[41]*draws[,41]
  ev <- out
  meanhelp <- means[43] + varmeans[1]*means[1]+ varmeans[2]*means[2]+varmeans[3]*means[3]+varmeans[4]*means[4]+varmeans[5]*means[5]+varmeans[6]*means[6]+varmeans[7]*means[7]+varmeans[8]*means[8]+adminvec[i]*means[9]+ (adminvec[i]^2)*means[10] + varmeans[41]*means[41]
  mean <- meanhelp
  lower <- quantile(ev,0.05)
  higher <- quantile(ev,0.95)
  meanvec <- c(meanvec,mean)
  lowervec <- c(lowervec,lower)
  highervec <- c(highervec,higher)
}

plots <- data.frame(adminvec,meanvec)
plotadmin <- data.frame(admin$adminpc[admin$adminpc<max])
colnames(plotadmin)[1] <- "adminpc"



#pdf("effect_square2.pdf",height=7, width=7*1.62)
ggplot(plots,aes(adminvec, meanvec))+stat_smooth(geom = "line")+
  theme_bw()+
  theme(axis.text=element_text(size=17),axis.title=element_text(size=17,face="bold"))+
  xlab("N. Top Tier Regional Govs, per 1m Citizens")+ylab("E(Expanded Services Index)")+geom_ribbon(aes(ymin=lowervec, ymax=highervec),alpha=0.2)+
  geom_point(data=plotadmin, aes(x=plotadmin$adminpc, y=min(lowervec)), alpha=0.2)+
  scale_x_continuous(breaks=c(0,1,2,3,4,5,6,7,8,9,10))+
  scale_y_continuous(breaks=c(-1,-0.5,0,0.5,1))+
  geom_vline(xintercept =0.28,colour="black",linetype = "solid")+
  annotate("text", label = "Nigeria", x =0.8, y = 0.75, size = 7, colour = "black")+
  geom_vline(xintercept =1.91,colour="black",linetype = "longdash")+
  annotate("text", label = "Uganda", x =2.5, y = 0.75, size = 7, colour = "black")+
  geom_vline(xintercept =2.32,colour="black",linetype = "dotted")+
  annotate("text", label = "Malawi", x =2.85, y = 0.5, size = 7, colour = "black")

################################################################
# Section F: National Leaders
################################################################

#####################
##### Figure F1 #####
#####################

PrezCandidates=read.dta("~PrezCandidates.dta", convert.factors=FALSE)
str(PrezCandidates)

Govs= PrezCandidates %>%
  select(year, Govenor) %>%
  filter(!is.na(Govenor), year>1962) %>%
  group_by(year) %>%
  summarise(x1 = round(mean(Govenor), digits=2))

#pdf(file="~ShareGovernors.pdf", width=10) 
ggplot(Govs, aes(x=year, y=x1)) + 
  geom_point(size=4, colour="red") + theme_bw() + geom_smooth(method = "loess", size = 1.5) +
  ggtitle("Governors' Share of Presidential Candidates (Africa)") + scale_x_continuous(name="") +
  scale_y_continuous(name="Share Governors") +
  theme(legend.position=c(0.8, 0.1),
        panel.grid.major.y = element_line(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"),
        strip.text = element_text(size=17),
        plot.title = element_text(size = 20, face = "bold", vjust = 1.5),
        axis.title.y = element_text(colour="grey20",size=15,face="bold"),
        axis.text.x = element_text(colour="grey20",size=15,face="bold"),
        axis.text.y = element_text(colour="grey20",size=15,face="bold"),  
        axis.title.x = element_text(colour="grey20",size=15,face="bold"), 
        axis.ticks.y = element_blank())

